### REPLICATION FILE -- APPENDIX, US
### Jonathan Homola & Margit Tavits
### "Contact Reduces Immigration-Related Fears for Leftist
### but Not for Rightist Voters"
### Comparative Political Studies

## clear environment, set seed, install/load packages
rm(list=ls())
set.seed(12435); options(stringsAsFactors=F)
# install.packages("stargazer"); install.packages("foreign")
# install.packages("psych"); install.packages("xtable")
library("foreign"); library("stargazer"); library("psych"); library("xtable")

## set working directory and source interaction functions
setwd(" ... ")
source("HomolaTavits_ContactFears_interactions.R")

## read in the main dataset 
data <- read.dta("HomolaTavits_ContactFears_US.dta")

## create contact factor score
contact <- prcomp(~ as.numeric(contact1) + as.numeric(contact2), 
                  data=data, na.action = na.exclude)
data$contact <- -contact$x[,1]

## Eigenvalue of first and second factor
(contact$sdev^2*2/sum(contact$sdev^2))[1]  ## 1.56
(contact$sdev^2*2/sum(contact$sdev^2))[2]  ## 0.44

## recode outcome variables
data$violence_com <- 5-as.numeric(data$violence_com)
data$violence_nat <- 5-as.numeric(data$violence_nat)
data$econ_hh <- 5-as.numeric(data$econ_hh)
data$econ_nat <- 5-as.numeric(data$econ_nat)
data$nat_id <- 5-as.numeric(data$nat_id)
data$amer_cult <- 5-as.numeric(data$amer_cult)

## create threat factor score
threat <- prcomp(~ violence_com + violence_nat + econ_hh + econ_nat + nat_id + amer_cult,
                 data=data, na.action = na.exclude)
data$threat <- threat$x[,1]

## Eigenvalue of first and second factor
(threat$sdev^2*6/sum(threat$sdev^2))[1]  ## 4.58
(threat$sdev^2*6/sum(threat$sdev^2))[2]  ## 0.56




#### Table OA3.1: Descriptive Statistics
data <- subset(data, !is.na(weights_sep))
data$female <- as.numeric(data$female)-1
data$midwest <- NA
data$midwest[!is.na(data$region)] <- 0
data$midwest[data$region==2] <- 1
data$south <- NA
data$south[!is.na(data$region)] <- 0
data$south[data$region==3] <- 1
data$west <- NA
data$west[!is.na(data$region)] <- 0
data$west[data$region==4] <- 1
data$Republican <- NA
data$Republican[!is.na(data$PID3)] <- 0
data$Republican[data$PID3==2] <- 1
data$Democrat <- NA
data$Democrat[!is.na(data$PID3)] <- 0
data$Democrat[data$PID3==1] <- 1
xtable(describe(data[,c("violence_com", "violence_nat",
                        "econ_hh", "econ_nat",
                        "nat_id", "amer_cult",
                        "contact", "threat",
                        "Republican", "Democrat", "PID7",
                        "female", "age", "education", "income",
                        "midwest", "south", "west")])[,c(2,3,4,8,9)])




#### Table OA3.2: Threat Factor Score
xtable(cor(cbind(data$violence_com, data$violence_nat, data$econ_hh,
                 data$econ_nat, data$nat_id, data$amer_cult, data$threat),
           use="complete.obs"))




#### Table OA3.3: Contact Factor Score
xtable(cor(cbind(as.numeric(data$contact1)*(-1), 
                 data$contact2, data$contact), use="complete.obs"))




#### Table OA4.1: The Effect of Contact on Immigration-Related Threat, the U.S.
mod0 <- lm(threat ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod1 <- lm(violence_com ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod2 <- lm(violence_nat ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod3 <- lm(econ_hh ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod4 <- lm(econ_nat ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod5 <- lm(nat_id ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod6 <- lm(amer_cult ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
stargazer(mod0, mod1, mod2, mod3, mod4, mod5, mod6,
          omit.stat=c("f", "ser", "bic", "ll", "adj.rsq"), 
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", font.size="small", digits=2, 
          dep.var.caption="Outcome variable:",
          dep.var.labels = c("Threat Factor", "Violence Com", "Violence Nat",
                             "Econ HH", "Econ Nat",
                             "Identity", "Culture"))




#### Table OA4.2: The Effect of Contact and Left-Right Affinity (binary) on
#### Immigration-Related Threat, the U.S.
data$PID3 <- relevel(as.factor(data$PID3), ref=3)
data2 <- subset(data, PID3!=3)
mod0 <- lm(threat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod1 <- lm(violence_com ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod2 <- lm(violence_nat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod3 <- lm(econ_hh ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod4 <- lm(econ_nat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod5 <- lm(nat_id ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
mod6 <- lm(amer_cult ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2)
stargazer(mod0, mod1, mod2, mod3, mod4, mod5, mod6,
          omit.stat=c("f", "ser", "bic", "ll", "adj.rsq"), 
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", font.size="small", digits=2, 
          dep.var.caption="Outcome variable:",
          dep.var.labels = c("Threat Factor", "Violence Com", "Violence Nat",
                             "Econ HH", "Econ Nat",
                             "Identity", "Culture"))




#### Table OA4.3: The Effect of Contact and Left-Right Affinity (continuous) on
#### Immigration-Related Threat, the U.S.
data$PID7 <- data$PID7-1
mod0 <- lm(threat ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod1 <- lm(violence_com ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod2 <- lm(violence_nat ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod3 <- lm(econ_hh ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod4 <- lm(econ_nat ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod5 <- lm(nat_id ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
mod6 <- lm(amer_cult ~ contact*PID7 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data)
stargazer(mod0, mod1, mod2, mod3, mod4, mod5, mod6,
          omit.stat=c("f", "ser", "bic", "ll", "adj.rsq"), 
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", font.size="small", digits=2, 
          dep.var.caption="Outcome variable:",
          dep.var.labels = c("Threat Score", "Violence Com", "Violence Nat",
                             "Econ HH", "Econ Nat",
                             "Identity", "Culture"))




#### Figure OA4.1: Marginal Effect of Contact on Threat Perceptions, the US
#pdf("FigOA4.1.pdf", height=8,width=7)
layout(matrix(c(0,0,0,1,1,1,1,1,0,0,0,
                2,2,2,2,2,0,3,3,3,3,3,
                4,4,4,4,4,0,5,5,5,5,5,
                6,6,6,6,6,0,7,7,7,7,7), 4, 11, byrow = TRUE))
interaction_plot_continuous(mod0, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.8,
                            title="Threat (score)")
interaction_plot_continuous(mod1, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="Violence, community")
interaction_plot_continuous(mod2, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="Violence, national")
interaction_plot_continuous(mod3, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="Economy, household")
interaction_plot_continuous(mod4, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="Economy, national")
interaction_plot_continuous(mod5, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="National identity")
interaction_plot_continuous(mod6, "contact", "PID7", "contact:PID7",
                            rugplot=F, histogram=F, pid=T,
                            ylabel="Marginal Effect of Contact",
                            xlabel="",
                            ymax=0.15, ymin=-0.4,
                            title="American culture")
par(mfrow=c(1,1))
#dev.off()




#### Table OA4.4: The Effect of Contact and Left-Right Affinity (binary) on
#### Immigration-Related Threat, the U.S. (Weighted)
moda <- lm(threat ~ contact + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data, weights=weights_sep)
mod0 <- lm(threat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod1 <- lm(violence_com ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod2 <- lm(violence_nat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod3 <- lm(econ_hh ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod4 <- lm(econ_nat ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod5 <- lm(nat_id ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
mod6 <- lm(amer_cult ~ contact*PID3 + female + age + as.numeric(education) + as.numeric(income) + as.factor(region), data2, weights=weights_sep)
stargazer(moda, mod0, mod1, mod2, mod3, mod4, mod5, mod6,
          omit.stat=c("f", "ser", "bic", "ll", "adj.rsq"), 
          star.cutoffs = c(0.05, 0.01, NA), 
          column.sep.width="1pt", font.size="small", digits=2, 
          dep.var.caption="Outcome variable:",
          dep.var.labels = c("Threat Factor", "Violence Com", "Violence Nat",
                             "Econ HH", "Econ Nat",
                             "Identity", "Culture"))